/**
 * @file   main8.cpp
 * @author Carlos Hernández García <carloshergar@usal.es>, Camilo Ruiz, Enrique Conejero
 * @date   Fri Apr  10 2009  
 *
 * Copyright (c) 2009 Carlos Hernández García
 *
 * @brief  Main program, in which we have implemented the moving ions
 *
 * 
 * Carlos Hernández García                          \n
 * <carloshergar@usal.es>                           \n
 * Universidad de Salamanca                         \n
 * Plaza de la Merced s/n                           \n
 * E-37008 Salamanca, SPAIN                         \n
 */
/********************************************************************************/


#include <iostream>
#include <vector>
#include <string>
//#include <cstring>

#include <alloc.h>
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <complex>
#include <fstream.h>

#include "Parameter.h"
#include "carpetas.h"
#include "gridobject.h"
//#include "../include/laserfield.h"
#include "routinesO.h"
#include "constants.h"

#include <fstream>


#define complex complex<double>

using namespace std;
using std::cout;
using std::vector;
using std::string;

//Parameter<string> jobname("jobname","Rana-pic: Main8");
//Parameter<double> lambda_nm("lambda_nm",32, "Wavelength in nm");
//                                          ^^ optional comment
//Parameter<vector<int> > point("point", vector<int>(3,0) );

Parameter<int>    np0("np0",4000,"Particle number"); 
Parameter<int>    ndg0("ndg0",8192,"Spatial grid points"); 
//Parameter<int>    ndt0("ndt0",32768,"Temporal grid points");
Parameter<int>    factorc_ndt0("factorc_ndt0",1,"Factor that determines ndt0 from ndg0 to verify deltag=c*deltat (factor=ndt0/ndg0) ");

//Mode of calculation
Parameter<int>    Fg_FieldsCalculation("Fg_FieldsCalculation",1,"Mode of calculation of the fields (1-Maxwell, 0-Potentials)");
Parameter<int>    Fg_ChargeCalculation("Fg_ChargeCalculation",0,"Mode of deposition of the charge (1-Lichters, 0-Enrique)");
Parameter<int>    Fg_CurrentCalculation("Fg_CurrentCalculation",0,"Mode of deposition of the current (1-Lichters, 0-Enrique)");
Parameter<int>    Fg_MovingIons("Fg_MovingIons",0,"If the ions are quit or not (1-movables, 0-quiet)");
Parameter<int>    Fg_InitialVelocity("Fg_InitialVelocity",0,"Initial Velocity (0-zero, 1-Birdshall, 2-Lichters)");


//Field parameters
Parameter<double> w0("w0",0.05,"Laser frequency");
Parameter<double> E0("E0",4,"Laser Amplitude");
Parameter<int>    nw("nw",20,"Number of laser cycles");
Parameter<int>    pulse_shape("pulse_shape",0,"Pulse shape (0-sin^2)");


//Temperature
Parameter<double> Etherm("Etherm",50/27.21,"Initial thermal energy of the particles");

//Density
Parameter<double> density("density",15.,"Density in critical density units");
Parameter<double> ionmass("ionmass",1386000.,"Mass of the ions");

//Thickness (in lambdas)
Parameter<double> thickness_lambda("thickness_lambda",1.,"Thickness of the material, in lambda units");

//Parameter for the average
Parameter<int> na("na",3,"Smoother");


//Output folders
//Parameter<string> foldername("foldername","pic");
Parameter<int> dir_limit("dir_limit",50,"Output parameter: Maximum number of directories that can be created");

//Output parameters
Parameter<int> field_seeing("field_seeing",4,"Output parameter: Where to see the field transmitted and reflected, in lambda units");
Parameter<int> sample_ndt("sample_ndt",20,"Output parameter: number of temporal samples per optical cicle");
Parameter<int> sample_ndg("sample_ndg",20,"Output parameter: number of spatial samples per wavelength");
Parameter<int> sampled_particles("sampled_particles",1000,"Output parameter: number of particles that are sampled");







int main(int argc, char *argv[])
{

  // read configuration file
  ConfigFileParser cfg( argv[1] );
    
  ParameterMap pMap;
  // initialize all Parameter<T> type objects from the config file
  pMap.init( cfg );
  
  
  //cout<<"Wavenumber: "<< 1/( lambda_nm * 1e-9 ) <<"\n";
  
  // string/map access:
  //double x; cfg.getValueInto(x, "lambda_nm"); // ConfigFileParser style
  //double y = pMap["lambda_nm"];               // ParameterMap style
  //string s = pMap["jobname"];                 // kind of "chameleon"
  
 /* try {
    cout << "Trying to pass a Parameter<double> to an int:\n";
    int i = pMap["lambda_nm"];
  } catch ( std::runtime_error e ) {
    cerr << "   ->    " << e.what() <<"\n";
  }*/
  
  cout<<pMap;
  
  /*************************************************/
  //Create folders and files.
  FILE *fp0; //Fichero de prueba
  FILE *fo1, *fo2, *fo3, *fo4; //Ficheros de salida de datos
  //string hh=pMap["foldername"];
  //char nomcarpeta=hh.c_str();
  
  //Creating a new data directory
  
  if(CrearCarpeta(dir_limit)==0)
    {	
      printf("Has superado el limite de directorios que tu pusiste, cazurro!!");
      exit(1);
    }
  
  
  fo1=fopen("output1.txt", "wt+"); //Fichero donde se almacenan los campos.
  if(fo1==NULL)
    {
      printf("Error opening the file.\n");
      exit(2);
    }
  fo2=fopen("output2.txt", "wt+"); //Fichero donde almaceno las cargas.
  if(fo2==NULL)
    {
      printf("Error opening the file.\n");
      exit(2);
    }
	fo3=fopen("output3.txt", "wt+"); //Fichero donde almaceno las densidades.
	if(fo3==NULL)
    {
      printf("Error opening the file.\n");
      exit(2);
    }
	fo4=fopen("output4.txt", "wt+"); //Fichero donde almaceno las densidades.
	if(fo4==NULL)
    {
      printf("Error opening the file.\n");
      exit(2);
    }
  
  
   ofstream myfile ("parameters.txt");
  if (myfile.is_open())
  {
      //  myfile << "This is a line.\n";
      //myfile << "This is another line.\n";
      myfile << pMap;
    //myfile.close();
  }
  else cout << "Unable to open file";
  
  
  /*************************************************/
  // Define variables
  /*************************************************/
  
  
  
	//int iaux,npion,npim;
	int npion=0;
  
	int np=pMap["np0"];//4000; //40000;//400;//3; //4000;              ///< Particle number
	int ndg=pMap["ndg0"];//1024*8; //1024*8; //200; //1024*8;           ///< Number of points in the spatial grid
	int factorc_ndt=pMap["factorc_ndt0"];//1024*32; //1024*32; //5000;//1024*32;      ///< Number of points in the temporal grid
	int ndt=ndg*factorc_ndt;
	
	/*************************************************/
	// Define transformation
	/*************************************************/
	double angle=45.;//_in=45.;
	
	double Beta        = sin( Pi/180 * angle );
	double Gamma       = 1.0 / cos( Pi/180 * angle );
	
	 cout << " Angle In " << angle << "\n";
	 cout << " beta " << Beta << "\n";
	 cout << " Gamma "<< Gamma << "\n";
	
	/*************************************************/
	// Define transformation
	/*************************************************/
	
	
	//Moving ions ATTENTION!! np must be even!!
	int npe; //Number of electrons and ions (the same)
	if(Fg_MovingIons==0) npe=np;
	else npe=np/2;
	
	myfile << "npe "<< npe<< "\n";
	
	//Field parameters
	laser laser0(pMap["w0"],pMap["nw"],pMap["E0"]);
	
	//Calculation parameters.
	//double lg=(laser0.nw+1)*laser0.lambda/4; //WATCH Fix here to make deltag=c*deltat
	//double tmax=(laser0.nw+1)*2*Pi/laser0.w0;
	double lg=(laser0.nw+5)*laser0.lambda/factorc_ndt; //WATCH Fix here to make deltag=c*deltat
	double tmax=(laser0.nw+5)*laser0.lambda/lightC_au;
	double tciclo=2*Pi/laser0.w0;
	double deltag=lg/ndg;
	double deltat=tmax/ndt;
	
	myfile << "ndt "<< ndt<< "\n";
	myfile << "lg "<< lg<< "\n";
	myfile << "tmax "<< tmax<< "\n";	
	myfile << "tciclo "<< tciclo<< "\n";	
	myfile << "deltag "<< deltag<< "\n";
	myfile << "deltat "<< deltat<< "\n";	
	
	//Material parameters
	double thickness_lambda_=pMap["thickness_lambda"];
	double thickness=thickness_lambda_*laser0.lambda;
	double xgleft=lg/2.-thickness/2.;
	double xgright=xgleft+thickness/2.;
	double nc=pow(Gamma,-2)*laser0.w0*laser0.w0*mass/4./Pi/q/q;  //Volumetric density

	double density_b=pow(Gamma,3) * density;//input.n_ion_over_nc / input.ppc[part->species];
	
	double n0=density_b*nc;//15*nc;
	
	
	double cells_per_wl=ndg/(lg/laser0.lambda); //cells per wavelength
	double vtherm=sqrt(density_b)*2*Pi/cells_per_wl*lightC_au;					
	
	//Parameter for the current calculation
	double current_factor=nc*thickness/deltag/npe; 
	double CalculaRofactor=n0*thickness/npe;

	//Limits of the material sheet
    int ileft=int(xgleft/deltag);
	int iright=int(xgright/deltag);

	myfile << "thickness "<< thickness<< "\n";
	myfile << "xgleft "<< xgleft<< "\t"<< " Left position of the sheet"<< "\n";	
	myfile << "xgright "<< xgright<< "\t"<< " Right position of the sheet"<< "\n";	
	myfile << "ileft "<< ileft<< "\t"<< " Left grid position of the sheet"<< "\n";	
	myfile << "iright "<< iright<< "\t"<< " Right grid position of the sheet"<< "\n";	
	myfile << "nc "<< nc<< "\n";
	myfile << "n0 "<< n0<< "\t"<< " (density*nc)"<< "\n";
	myfile << "lambda "<< laser0.lambda<< "\n";

	//Parameters for the output
	double sample_inc_t=tciclo/sample_ndt;
	double sample_inc_g=laser0.lambda/sample_ndg;
	double sample_inc_part=np/sampled_particles;
	
	myfile << "sample_inc_t "<< sample_inc_t<< "\n";
	myfile << "sample_inc_g "<< sample_inc_g<< "\n";
	myfile << "sample_inc_part "<< sample_inc_part<< "\n";
	myfile << "cells_per_wl "<< cells_per_wl<< "\n";
	myfile << "vtherm "<< vtherm<< "\n";

	
				
	particle *part;
	part=new particle[np];
	cell1d *cell;
	cell=new cell1d[ndg];
	
		
	double xg;//,deltag;//deltat,tmax,gamma;
	double Epx,Epy,Epz,Bpx,Bpy,Bpz;
	
	//FFT variables
	double sp1[ndt],sp2[ndt],sp3[ndt],sp4[ndt];//,raux;
	 
	//Particle parameters
	//vector<double> qp(np,0.),mp(np,0.);
	
	int iflagion=0;
	int labelx=1, labely=2, labelz=3; //labels used as arguments for the function CalculaJ
	double igamma; //just for ptinting in output2.txt
	//Particles that are rejected in the other side
	int coolparticles=0;
	
	printf("deltat=%f, deltag=%f, lg=%f, tmax=%f, lambda=%f, deltag/deltat=%f\n", deltat, deltag, lg, tmax, laser0.lambda, deltag/deltat);
		

	/******************************************************/
	/* Initialize variables
	/* Calculate Densities, Fields and position at the beggining
	/******************************************************/
	

	if (iflagion==0)
	{
		npion=np;
	}
	else
	{
		npion=0;
	}
	
	
	//Initialization of the position of the cells
	for(int l=0; l<ndg; l++)
	{	
		cell[l].vx.x=(l)*deltag;	
		cell[l].ct.x=(l)*deltag+(deltag/2);
	}

	//Initialization of the particles	
	if(iflagion==0)
	{
		InicioVelocidades(part, Fg_InitialVelocity, Etherm, npion, vtherm, fo2);

		for(int i=0;i<npe;i++)
		{
			//Electrons
			part[i].zn=density_b;
			part[i].xp=thickness/npe*1.*(i+1)+xgleft;
			part[i].xpold=thickness/npe*1.*(i+1)+xgleft;
			
			//Transforming velocities to the Boosted frame
			
			part[i].vpx    = part[i].vpx * sqrt(1.0-Beta*Beta) / ( 1 - part[i].vpy *Beta/lightC_au );
			part[i].vpz    = part[i].vpz * sqrt(1.0-Beta*Beta) / ( 1 - part[i].vpy *Beta/lightC_au);
			part[i].vpy    = ( part[i].vpy - Beta*lightC_au ) /  ( 1 - part[i].vpy *Beta/lightC_au);


			//Converting velocities into momentums	
			double gamma=1./sqrt(1.-(part[i].vpx*part[i].vpx+part[i].vpy*part[i].vpy+part[i].vpz*part[i].vpz)/lightC_au/lightC_au);
			part[i].vpx=part[i].vpx*gamma;
			part[i].vpy=part[i].vpy*gamma;
			part[i].vpz=part[i].vpz*gamma;
			
			if(Fg_MovingIons)//Ions
			{	part[i+npe].zn=density_b;
				part[i+npe].q=1.;
				part[i+npe].m=ionmass;
				part[i+npe].qm=part[i+npe].q/part[i+npe].m;
				part[i+npe].xp=thickness/npe*1.*(i+1)+xgleft;
				part[i+npe].xpold=thickness/npe*1.*(i+1)+xgleft;

				//Transforming velocities to the Boosted frame
				part[i+npe].vpx    = part[i+npe].vpx * sqrt(1.0-Beta*Beta) / ( 1 - part[i+npe].vpy *Beta/lightC_au );
				part[i+npe].vpz    = part[i+npe].vpz * sqrt(1.0-Beta*Beta) / ( 1 - part[i+npe].vpy *Beta/lightC_au);
				part[i+npe].vpy    = ( part[i+npe].vpy - Beta*lightC_au ) /  ( 1 - part[i+npe].vpy *Beta/lightC_au);
								
				gamma=1./sqrt(1.-(part[i+npe].vpx*part[i+npe].vpx+part[i+npe].vpy*part[i+npe].vpy+part[i+npe].vpz*part[i+npe].vpz)/lightC_au/lightC_au);
				part[i+npe].vpx=part[i+npe].vpx*gamma;
				part[i+npe].vpy=part[i+npe].vpy*gamma;
				part[i+npe].vpz=part[i+npe].vpz*gamma;
			}
		}
	}//end if
	
	//Printing the initial distribution of velocities
	int l;
	for(l=0; l<npion; l++)
	{
		igamma=part[l].igamma();
		fprintf(fo2, "%e\t%e\t%e\t%e\n", 0*deltat/tciclo, part[l].vpx*igamma, part[l].vpy*igamma, part[l].vpz*igamma);
	}
	
	//CHarge and Current deposition
	if(Fg_ChargeCalculation==0)
		CalculaRo(part,cell,CalculaRofactor, deltag, na, np, ndg, fp0);
	else
		for(int k=0; k<np; k++)
			Charge_deposition (cell, part[k], deltag, current_factor);
	
	if(Fg_CurrentCalculation==0)
	{
		CalculaJ(labelx,part,cell,CalculaRofactor, deltat, deltag, na, np, ndg, fp0, fp0); //labelx is 1 and means 'x'
		CalculaJ(labely,part,cell,CalculaRofactor, deltat, deltag, na, np, ndg, fp0, fp0); //labely is 2 and means 'y'
		CalculaJ(labelz,part,cell,CalculaRofactor, deltat, deltag, na, np, ndg, fp0, fp0); //labelz is 3 and means 'z'
	}
	else
	{	printf("Calculo como Lichters la corriente. npion=%d\n", npion);
		for(int k=0; k<npion; k++)
		{	//printf("Llamada a la funcion corriente\n");
			Current_deposition (cell, part[k], deltag, current_factor, 1, fo2, ndg);
		}
	}		
			
	
	
	
	/******************************************************/
	/* Central part of the program
	/* Time and particle loop 
	/******************************************************/
	
	//Parameters for sampling the output
	double simulationtime=0;
	double simulationspace=0;
	double simulationpart=0;
	double simulationsteps_g=0;
	double simulationsteps_t=0;
	double simulationsteps_p=0;
	
	for(int n=0;n<ndt;n++)
	{
		if(n%100==0) 
		{
			printf("n=%d\n", n);
		}
	
		double t=n*deltat;
		simulationtime+=deltat;
		
		for(int j=0;j<npion;j++)
		{
			interpola(part[j].xp,cell,Epx,Epy,Epz,Bpx,Bpy,Bpz,t,ndg,deltag,na, fp0);
		
			//Boris algorithm to solve the equation of motion
			ParticleMover(part[j],Epx,Epy,Epz,Bpx,Bpy,Bpz,deltat,lg, coolparticles, n, fo2);
		
		}//End of the Boris algorithm to solve the equation of motion
		if(Fg_ChargeCalculation==0)
		{
			CalculaRo(part,cell,CalculaRofactor, deltag, na, np, ndg, fp0);
		}		
		else
		{	
			for (int l=0;l<ndg;l++)
				cell[l].ct.Rox=0.;
	
			for(int k=0; k<np; k++)
			{	
				Charge_deposition (cell, part[k], deltag, current_factor);
			}
		}
		if(Fg_CurrentCalculation==0)
		{
			CalculaJ(labelx,part,cell,CalculaRofactor, deltat, deltag, na, np, ndg, fp0, fp0); //labelx is 1 and means 'x'
			CalculaJ(labely,part,cell,CalculaRofactor, deltat, deltag, na, np, ndg, fp0, fp0); //labely is 2 and means 'y'
			CalculaJ(labelz,part,cell,CalculaRofactor, deltat, deltag, na, np, ndg, fp0, fp0); //labelz is 3 and means 'z'
		}
		else
		{
			for (int l=0;l<ndg;l++)
			{
				cell[l].vx.Jx=0.;
				cell[l].ct.Jy=0.;
				cell[l].ct.Jz=0.;
			}
			for(int k=0; k<npion; k++)
			{	
				Current_deposition (cell, part[k], deltag, current_factor, n, fo2, ndg);
			}
			
		}
		//printf("Rox[205]=%f\tRoy[205]=%f\tRoz[205]=%f\tRox[298]=%f\tRoy[298]=%f\tRoz[298]=%f\n", cell[297].ct.Rox, cell[297].ct.Roy, cell[297].ct.Roz, cell[298].ct.Rox, cell[298].ct.Roy, cell[298].ct.Roz);
		//printf("%d, Rox[297]=%e\tJx[297]=%e\tJy[297]=%e\tJz[297]=%e\tRox[298]=%e\tJx[298]=%e\tJy[298]=%e\tJz[298]=%e\n",n, cell[297].ct.Rox, cell[297].vx.Jx, cell[297].ct.Jy, cell[297].ct.Jz,cell[298].ct.Rox, cell[298].vx.Jx, cell[298].ct.Jy, cell[298].ct.Jz);
		if(n==1)
		{	float Q=0.;
			for(int i=0; i<ndg; i++)
			{	//printf("[%d]=%e\n",i, cell[i].ct.Rox);
				Q+=cell[i].ct.Rox;
			}
			printf("QTOTAL=%e\n", Q);
		}
		
		
                     
		/******************************************************/
		/* Calculate Ex, Ey, Ez, Bx, By, Bz, directly
		/* By solving the Maxwell equations 
    	/******************************************************/
           
        CalculaCamposInternos(cell, deltat, xgleft, laser0, ndg, deltag, fo1);
       	for(int k=0; k<ndg; k++)
        {   cell[k].vx.Ey+=E0y(k*deltag,t,xgleft,laser0, pulse_shape);
				
            cell[k].vx.Bz+=B0z(k*deltag,t,xgleft,laser0, pulse_shape);
        }
		//printf("\tE0y[205]=%e\n", E0y(205*deltag,t,xgleft,laser0));
		//printf("Ex[297]=%f\tEy[297]=%f\tEz[297]=%f\tEx[298]=%f\tEy[298]=%f\tEz[298]=%f\n", cell[297].vx.Ex, cell[297].vx.Ey, cell[297].vx.Ez, cell[298].vx.Ex, cell[298].vx.Ey, cell[298].vx.Ez);
	
		
		if(iflagion!=0)
		{
			//Ionization();
			//Calculap0();
		}
		
		
		//Calculate the transmitted and reflected field.
		sp1[n]=cell[int(((xgleft-field_seeing*laser0.lambda)/deltag)-1)].vx.Ey-E0y( (xgleft-field_seeing*laser0.lambda)-1,t+deltat,xgleft,laser0, pulse_shape);//wp,kp,E0,k0,w0);
		sp2[n]=cell[int(((xgright+field_seeing*laser0.lambda)/deltag)-1)].vx.Ey;
		
		
		
		
		
		
		/**********************************WRITING THE OUTPUT***************************/
		
		fprintf(fo1, "%d\t%e\t%e\t%e\t%e\t%e\t%e\t%e\t%e\n", n, t/2/Pi*laser0.w0, sp1[n], sp2[n],E0y(xgleft-field_seeing*laser0.lambda,t+deltat,xgleft,laser0,pulse_shape), cell[ileft].vx.Ex, cell[iright].vx.Ex, cell[ileft].ct.Rox, cell[iright].ct.Rox);
			
		//if(n%(sample_division_ndt)==0)
		simulationspace=0;
		simulationsteps_g=0;
		if(simulationtime>=sample_inc_t)
		{	
		
			//fprintf(fo1, "%d\t%e\t%e\t%e\t%e\t%e\t%e\t%e\t%e\n", n, t/2/Pi*laser0.w0, sp1[n], sp2[n],E0y(xgleft-field_seeing*laser0.lambda,t+deltat,xgleft,laser0,pulse_shape), cell[ileft].vx.Ex, cell[iright].vx.Ex, cell[ileft].ct.Rox, cell[iright].ct.Rox);
			
			simulationpart=0;
			simulationsteps_p=0;
			for(int k=0; k<npion; k++)
			{   
				if (npion>sampled_particles)
				{
					simulationpart++;
					if(simulationpart>=sample_inc_part)
					{
						igamma=part[k].igamma();
						fprintf(fo2, "%e\t%e\t%e\t%e\n", n*deltat/tciclo, part[k].vpx*igamma, part[k].vpy*igamma, part[k].vpz*igamma);
						simulationpart=simulationpart-sample_inc_part;
						simulationsteps_p++;
					}
				}
				else
				{	
					igamma=part[k].igamma();
					fprintf(fo2, "%e\t%e\t%e\t%e\n", n*deltat/tciclo, part[k].vpx*igamma, part[k].vpy*igamma, part[k].vpz*igamma);
					simulationsteps_p++;
				}
			}
			for(int k=0; k<ndg; k++)
			{   
				simulationspace+=deltag;
				if (simulationspace>=sample_inc_g)
				{
					fprintf(fo3, "%e\t%e\t%e\n", k*deltag/laser0.lambda, n*deltat/tciclo, cell[k].ct.Rox/nc);
					fprintf(fo4, "%e\t%e\t%e\n", k*deltag/laser0.lambda, n*deltat/tciclo, cell[k].vx.Ex);
					simulationspace=simulationspace-sample_inc_g;
					simulationsteps_g++;
				}
			}
			simulationtime=simulationtime-sample_inc_t;
			simulationsteps_t++;
		}
				
		
	}//end time loop
	 
	printf("Coolparticles=%d\n", coolparticles);
	
	myfile << "Coolparticles "<< coolparticles<< "\n";
	myfile << "simulation_steps_t "<< simulationsteps_t<< "\n";
	myfile << "simulation_steps_g "<< simulationsteps_g<< "\n";
	myfile << "simulation_steps_p "<< simulationsteps_p<< "\n";


	myfile.close();
	
	fclose(fo1);
    fclose(fo2);
	fclose(fo3);
	fclose(fo4);


	
	
}

void f1(grid a)
{
	cout << a.dx1;
}
